Exercise 12.1

The goal here is to create a prediction model for the value of benthic over the Chesepeake bay region.

This prediction model will be graphically portrayed with a:

1. Heatmap
   1. with base plot
   2. wigh ggplot2 over an actual map
2. Contour plot
   1. with base plot
   2. With ggplot2 over an actual map
3. 3d Surface plot
   1. with base plot
   2. with a dynamic plotly surface plot

Preamble

Import Library

The following code will automatically install packages that are not already installed and load them into the library:

if(require('pacman')){
  library('pacman')
}else{
  install.packages('pacman')
  library('pacman')
}
pacman::p_load(sp, EnvStats, gstat, ggplot2, rmarkdown, reshape2, ggmap, RColorBrewer, parallel, dplyr, plotly)

Create Assignments

Note that Benthic.df is included in the EnvStats package and is not a data frame but an sp class, which is very similar to a data frame

#Assignments
benthic <- Benthic.df
lat   <- benthic$Latitude
lon   <- benthic$Longitude
index <- benthic$Index

Inspect the Data Set

head(benthic, 3)

Lattitude and Longitude

• Lattitude is a line that runs east to west that measures how far north of the equatour a location is
  • Lattitude is assigned to the y-axis
• Longitude is a line that runs north to south that measures east/west
  • Longitude is assigned to the x-axis

Assign the Co-ordinates

The Benthic Data set is an sp class objects and requires for its co-ordinates to be set:

coordinates(benthic) = ~Longitude + Latitude

Preliminary Step: Create a krige prediction Model

In order to create a krige prediction model two things are required:

• Fitted Variogram Model (e.g. an exponential model)
• A domain values over which to predict values

Create a variogram Model

Create a Variogram

If a 4th degree polynomial model was used to model the benthic values over space, a variogram can be made from that model thusly:

(The quadratic model is a required parameter of the kriging process, in exercise 12.2 a quadratic model will be used, the spacific mathematics of this though is somewhat beyond me)

x <- lon
y <- lat
benthic.vg <- variogram(object = index ~ 
                          x + y + 
                          x^2 + x*y + y^2 + 
                          x^3 + x^2*y + x*y^2 + y^3 +
                          x^4 + x^3*y + x^2*y^2 + x*y^3 + y^4,
                        data = benthic)
head(benthic.vg)

    np       dist     gamma dir.hor dir.ver   id
1 4545 0.02410565 0.8506164       0       0 var1
2 4149 0.06703996 0.9393502       0       0 var1
3 3883 0.11396396 1.1513831       0       0 var1
4 4023 0.16138943 1.4599238       0       0 var1
5 5963 0.20615519 1.5641015       0       0 var1
6 6732 0.24983031 1.8992994       0       0 var1

Fit an exponential model to that variogram

benthic.vg.fit.exp <- fit.variogram(benthic.vg, model=vgm(1,"Exp", 0.5,1))

Thus the variogram model needed for the krige prediction has been ascertained

Create a Grid of Domain values to use as predictor values

Create a sequence of incrementally increasing for the X and Y axis

This will essentially act as the resolution of the overlay, Most monitors are about 140 ppi, hence assuming the plot will be a 9 inch square, about 800 x 800 pixels should look right, but 200 x 200 will be the used due to performance.

xgrid       <-  seq(min(lon), max(lon), length.out = 200)
ygrid       <-  seq(min(lat), max(lat), length.out = 200)

Combine into a domain base

This domain base however has to be a spatial object, this can be acheived by defining co-ordinates in the data frame

xy.surface  <- expand.grid(lon = xgrid, lat = ygrid)
coordinates(xy.surface) = ~lon + lat
head(xy.surface, 3)

class       : SpatialPoints 
features    : 3 
extent      : -77.3157, -75.9, 38.0018, 39.4883  (xmin, xmax, ymin, ymax)
coord. ref. : NA 

Create the Krige prediction object

benthic.pred <- krige(formula = Index ~1,
                           locations =  benthic,
                           newdata = xy.surface,
                           model= benthic.vg.fit.exp)

[using ordinary kriging]

head(benthic.pred)

class       : SpatialPointsDataFrame 
features    : 6 
extent      : -77.3157, -77.28013, 38.0018, 38.0018  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 2
names       :        var1.pred,         var1.var 
min values  : 3.12118704456262, 2.29623742979447 
max values  : 3.14485115010402, 2.32162957155164 

Now there is a prediction model, if that prediction model is extrapolated over our domain (xy.surface) the various plots can be created

Part 1: Create a heatmap

In order to create a heatmap, it is first necessary to take the sp prediction object and make it into pixel data, otherwise the heatmap won’t look quite right:

class(benthic.pred)

[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

gridded(benthic.pred) = TRUE
class(benthic.pred)

[1] "SpatialPixelsDataFrame"
attr(,"package")
[1] "sp"

spplot(benthic.pred["var1.pred"],main="Predictions of Benthic Index - ordinary kriging")

Create it in ggplot2

First it is necessary to create a data frame from the sp object:

benthic.pred.df      <- cbind(benthic.pred@coords, benthic.pred@data)
benthic.pred.df.tidy <- melt(benthic.pred.df)

No id variables; using all as measure variables

head(benthic.pred.df)

ggplot() +
  geom_tile(data = benthic.pred.df, aes(lon, lat, fill = var1.pred)) +
  labs(fill = "Benthic Index") # 

This plot would be easier to read if it had discrete intervals rather than a continuous scale:

benthic.pred.df$bin <- cut(benthic.pred.df$var1.pred,
                            breaks = as.vector(seq(from = 0, to = 5, length.out = 20)))
ggplot() +
  geom_tile(data = benthic.pred.df, aes(lon, lat, fill = bin)) +
  labs(fill = "Benthic Index", 
       title = "Heatmap of Benthic Index Values",
       x = "Longitude",
       y = "Lattitude")

NA

Map overlay

Now by adjusting the colour scale and opacity, this can be used as an overlay for the map from exercise 10.1:

Part 2: Contour Map

In order to view the surface plot, contours can also be used.

Recall from Week 11/Exercise 10, that contour and surface plots work a little differently to ordinary plots.

A domain grid of equally spaced \(x\) and \(y\) values is required and a matrix of \(z\) values that can be overlayed onto that domain grid in order to provide the surface values required to create the surface plot.

The method from Lecture 10/Wk. 11 involved predicting over the square domain of xy.surface, this won’t work here because instead of the predict function, the krige function was used, hence it is necessary to convert the output from the krige function, into the output that would have been given by the predict function, which is the required \(z\) surface matrix.

This can be acheived with the data.matrix function:

z_surface_matrix <- data.matrix(benthic.pred[1])

Base Package

Using base R packages, a contour plot can be created with the contour function:

contour(xgrid, ygrid, z_surface_matrix,
          xlab="-Longitude (degrees West)", ylab="Latitude (degrees North)",
          main = "contour plot of Krige Prediction",
          labcex = 0.8)
  mtext("Krige using a 4th degree polynomial model", cex = 0.5)

Using ggplot2

In order to create a contour plot in ggplot 2, it will be necessary to use ‘tidy’ data, which is a table of data such that each column corresponds to a variable:

rownames(z_surface_matrix) <- xgrid
  colnames(z_surface_matrix) <- ygrid
  benthic.pred.melt <- melt(z_surface_matrix,
                            varnames = c("Longitude", "Latitude"),
                            value.name = "Benthic_Index")
  head(benthic.pred.melt)

Observe that this is another way of getting the data frame that was required in part 1 (benthic.pred.df), it is included here as another method and for the sake of completion, also it matches the method used in week 10/lecture 11.

Now using ggplot2 a contour plot, can be created:

ggplot(data = benthic.pred.melt, aes(x = Longitude,
                                       y = Latitude,
                                       z = Benthic_Index,
                                       colour = ..level..),
         size = 0.5) +
    labs(x = "Longitude", y = 'Latitude',
         title = 'Benthic Index',
         subtitle = "Chesapeake Bay, Maryland",
         caption = "Higher index value is better",
         col = 'Benthic Index', size = 'Benthic Index') + 
    geom_contour(lwd = 1) +
  theme_classic()

Map Overlay

The advantage to using ggplot2, is the map overlay can be taken advantage of:

Map Overlay with Heatmap

This can be combined with the heatmap layer:

  ggmap(ggmap = map) +
    labs(fill = "Benthic Index",
         x = "Longitude",
         y = "Lattitude", 
         title = "Heatmap of Benthic Index") + 
    guides(alpha=FALSE) +
    geom_tile(data = benthic.pred.df, aes(lon, lat, fill = bin, alpha = var1.pred)) +
    geom_contour(data = benthic.pred.melt, aes(x = Longitude,
                                               y = Latitude,
                                               z = Benthic_Index),
                col = 'indianred',
                binwidth = 0.2,
                size = 0.6) +
    scale_fill_manual(values = cols) 

Part 3: Surface Plots

A surface plot allows a three dimensional view of the model

Base

In base package, a surface plot can be made:

persp(xgrid, ygrid, z_surface_matrix,
        xlim = c(-77.3, -75.9),   #the values domain needs to be adjusted above
        ylim = c(38.1, 39.5),          #Make sure to set an appropriate domain
        # zlim = c(0, 6),
        theta = -45, phi = 30, d = 0.5,
        xlab="-Longitude (degrees West)",
        ylab="Latitude (degrees North)",
        zlab="Benthic Index", ticktype = "detailed") %>% invisible()

surface extends beyond the box

  title(main=paste("Surface Plot of Benthic Index", "Loess Smoothing", sep="\n"))   

Plotly

Plotly has the advantage of being fully interactive and not as difficult to constrain within a plot/box area:

p <- plot_ly(x = ygrid, y = xgrid, z = z_surface_matrix) %>%
    add_surface() %>%
    layout(
      title = "Benthic Index in  Chesapeake Bay",
      scene = list(
        xaxis = list(title = "Longitude",
                     range = c(38.2, 38.9)),  #You shouldn't need to edit the
        yaxis = list(title = "Latitude",      #The domain here, do it above
                     range = c(-77, -75.9)),
        zaxis = list(title = "Benthic Index")
      ))
                                                                  
  p

The surface plot really reinforces the fact that the bay has the lowest (i.e. best) benthic index value.